% This function calculates the nonlinear constraint as in Equation (D.38) when EIS = 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [c, ceq] = projection_constraint_PD_EIS_1(beta_PD_proj, I , v, fspace_tmudc, tmu_dc_grid)

% INPUTS: 
% -- beta_PD_proj: Coefficients of polynomials from P/D approximation
% -- I: Parameter structure for other parameters
% -- v: Number of points for intergral
% -- fspace_tmudc: Function space that contains all the projection
%                  information for d_t - c_t and \tmu_t
% -- tmu_dc_grid: Two-dimensional grid points to be evaluated for tmu_t and d_t - c_t
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

coeff_1 = sqrt(1 + I.nu) * I.sigma;
coeff_2 = coeff_1 * I.nu;
coeff_3 = I.lambda * coeff_1 - I.xi * I.sigma;

[epsilon_tmu,weight_tmu] = qnwlege(v,-8,8);             %For epsilon from -8 to 8, do Gaussian-Quad sampling
[epsilon_sigmad,weight_sigmad] = qnwlege(v,-8,8);       %For epsilon from -8 to 8, do Gaussian-Quad sampling


tmu_grid = cell2mat(tmu_dc_grid(1,1));   %Get the grid for tmu to evaluate at time t
dc_grid = cell2mat(tmu_dc_grid(1,2)); %Get the grid for d_t-c_t to evaluate at time t

% Degrees of approximation along each dimension
n(1) = size(tmu_grid,1);
n(2) = size(dc_grid,1);

H_matrix_new = zeros(n(1)*n(2), v^2);%Value of log P/D ratio at time t+1. x for grid pairs, y for shocks
for j = 1 : n(2)
    for i = 1 : n(1)
        tmu_grid_new_temp = tmu_grid(i) + coeff_2 * epsilon_tmu;
        dc_grid_new_temp = (1 - I.alpha) * dc_grid(j) + I.alpha * I.mu_dc + (I.lambda - 1) * (tmu_grid(i) + coeff_1 * epsilon_tmu) + I.sigma_d * epsilon_sigmad;
        grid_new_temp{1,1} = tmu_grid_new_temp;
        grid_new_temp{1,2} = dc_grid_new_temp;
        H_matrix_new(n(1) * (j - 1) + i , :) = (funbas(fspace_tmudc, grid_new_temp) * beta_PD_proj)'; 
    end
end

expH_matrix_new = exp(H_matrix_new);

%Functions related to log P/D ratio in the denominator
fbasis_dc = funbas(fspace_tmudc,tmu_dc_grid);
H_matrix = bsxfun(@times , fbasis_dc * beta_PD_proj , ones(1 , v^2));
expH_matrix = exp(H_matrix);

%Remaining shocks from \Delta d_{t+1} and SDF
Add_Shock_Matrix = exp(bsxfun(@times , ...
                       kron(ones(v,1) , coeff_3*epsilon_tmu) + kron(I.sigma_d * epsilon_sigmad , ones(v,1)) ,...
                       ones(1 , n(1) * n(2)) ) )';

%Shocks density function
Shock_density_matrix = (0.5/pi) * exp( -0.5 * ...
bsxfun(@times , kron(ones(v , 1) , epsilon_tmu).^2 + kron(epsilon_sigmad , ones(v,1)).^2 , ones(1 ,n(1) * n(2))) )';

%Calculate the expectation term as integral 
Expectation_integral = (Shock_density_matrix .* Add_Shock_Matrix .*((expH_matrix_new + 1)./expH_matrix)) * kron(weight_sigmad , weight_tmu);

%Multiply the expectation term by other time t terms
Multiplier_matrix =  exp(I.tmu_m + (I.lambda-1) * kron(ones(n(2) , 1) ,tmu_grid ) - I.alpha * kron(dc_grid,ones(n(1),1)) + I.alpha * I.mu_dc);

c = [];
ceq = Multiplier_matrix .* Expectation_integral - 1;

end